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ABSTRACT 

The  illumination  condition  of  an  artificial  earth  satellite 
is  investigated  in  relation  to  the  geographic  location  of  the 
sub- satellite  point.     A  stochastic  problem  is  formulated  by 
postulating  a  "random"  satellite,  and  the  solution  to  this 
problem  is  indicated  analytically.     A  statistical  solution 
obtained  from  a  digital  computer  simulation  program,  together 
with  a  sensitivity  analysis  of  parametric  variations. 


li 


TABLE  OF  CONTENTS 

Page 

List  of  Illustrations  iv 

Glossary  vii 

Chapter  I  -  Orbital  Mechanics  5 

Chapter  II  -  Probability  Formulation  21 

Chapter  III  -  Simulation  Problem  27 

Bibliography  30 

Appendix  A  -  Design  of  a  Computer  Simulation  Program  A-l 

Appendix  B  -  Summary  of  Program  Output  B-l 


in 


LIST  OF  ILLUSTRATIONS 
Figure  Title  Page 

1-1  -  Representation  of  Earth -Sun  vector  E                                          7 

1-2  -  Representation  of  Earth- Satellite  vector  r                              9 

1-3  -  Simplified  view  of  Earth-Satellite  vector  r                             10 

1-4  -  Location  of  sub- satellite  point  at  time  t                                  12 

1-5  -  Plane  containing  E  and  r  vectors                                                    14 

1-6  -  Mercator  projection  of  track,  first  case                                 17 

1-7  -  Mercator  projection  of  track,   second  case                            18 

A-l  -  Mercator  projection  of  track,  first  case  (modified)       A-2 

A-2  -  Mercator  projection  of  track,   second  case  (modified)  A-3 

A-3  -  Mercator  projection  of  track,   special  case                         A-6 

B-l  -  Prob.   of  passage  and  illumination,   L  =  25° 

B-2  -  Prob.   of  passage  and  illumination,   L  =  35° 

B-3  -  Prob.   of  passage  and  illumination,   L  -  45 

B-4  -  Prob.   of  passage  and  illumination,   L  =  55° 

B-5  -  Prob.   of  passage  and  illumination,  1-28    ,   L  =  25 

B-6  -  Prob.   of  passage  and  illumination,  I  -  36    ,   L  -  35 

o  o 

B-7  -  Prob.   of  passage  and  illumination,  I  =  44    ,   L  =  45 

B-8  -  Prob.   of  passage  and  illumination,  I  -  52    ,   L  =  55 

B-9  -  Prob.   of  passage  and  illumination,  I  -  60    ,   L  =  55 

B-10  -  Prob.   of  passage  and  illumination,  I  =  68°,   L  =  65 

B-ll  -  Prob.   of  passage,   L  =  25° 


IV 


B-12    -         Prob.   of  passage,   L.  =  35 


O 

B-13  -  Prob.  of  passage,  L  =  45 
B-14  -  Prob.  of  passage,  L  =  55 
B-  15  -    Prob.  of  passage  ,  L,  =  65° 


PREFACE 

The  general  purpose  of  this  thesis  is  to  investigate  certain  aspects 
of  the  illumination  of  artificial  earth  satellites  by  the  sun,  in  relation 
to  the  track  of  the  sub- satellite  point  as  measured  in  a  geographic 
coordinate  system.     Imposed  upon  this  study  is  a  stochastic  model 
of  the  earth- satellite  system. 

An  analytic  approach  to  the  basic  problem  was  made  ,  followed 
by  a  computer  simulation  program  for  a  statistical  solution.     The 
basic  problem  is  set  forth  in  the  Introduction. 

I  wish  to  express  my  appreciation  to  the  following  persons  for 
their  contributions  to  this  study: 

Mr.   E.   J.    Leinfelder,  Chance-Vought  Corp.  ,  Astronautics 
Division,  who  first  posed  the  problem  to  me,  and  who  provided  the 
stimulus  and  inspiration  leading  to  this  study; 

Professor  W.    P.   Cunningham,   U.   S.   Naval  Postgraduate  School, 
who  indicated  the  method  of  formulating  the  problem  in  a  manner 
which  led  to  the  desired  results; 

Professor  R.    M.    Thatcher  and  Miss  Pat  Hoang  ,   U.   S.   Naval 
Postgraduate  School,  who  assisted  me  in  no  small  measure  in  the  design, 
writing  and  de-bugging  of  the  computer  simulation  program;    and 

Professor  R.   R.   Read,   U.   S.   Naval  Postgraduate  School,  who 
detected  initial  errors  in  the  chapter  on  probability  formulation. 
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GLOSSARY 

Symbol  Definition 

a  semi-major  axis  of  orbital  ellipse 

dl  apparent  change  in  longitude  of  sub- satellite 

point  from  initial  reference  time  to  time  t 

h  altitude  of  satellite  above  surface  of  earth 

i  inclination  of  satellite  orbital  plane  to 

equatorial  plane 

i  '  inclination  of  equatorial  plane  to  ecliptic  plane 

W  longitude  of  designated  geographical  points 

1  '.  longitude  of  sub- satellite  point  at  time  t 

1'  longitude  of  ascending  node 

r  unit  vector  from  center  of  earth  to  satellite 

r1  projection  of  r    onto  equatorial  plane 

s  specified  angular  distance  of  sub- satellite  point 

from  designated  geographical  points 

t  instantaneous  time  after  initial  reference  time 

w  angular  rotational  velocity  of  earth  about  polar 

axis  in  radians  per  unit  time 

w  angular  velocity  of  satellite  in  radians  per  unit 

time 

A  angle  between  E     and  r    vectors 

C  angle  at  satellite  subtended  by  the  center  of  the 

earth  and  the  surface  of  the  earth 

D  angle  between  vernal  equinox  and  ascending  node 

of  satellite 


vn 


Symbol  Definition 

E  unit  vector  from  center  of  earth  to  sun 

F  angle  between  ascending  node  and  projection  of  r 
onto  equatorial  plane 

G  angle  between  ascending  node  and  r    vector 

G'  universal  gravitational  constant 

H  angle  between  vernal  equinox  and  E    vector 

L  latitude  of  designated  geographical  points 

L'  angle  between  r    and  r1 

M  event  that  desired  passage  occurs  during  first  orbit 

M'  mass  of  the  earth 

N  event  that  desired  illumination  occurs  at  time  of 
closest  approach  to  designated  geographical  point 

R  radius  of  the  earth 

T  orbital  period  of  satellite 
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INTRODUCTION 

It  is  proposed  to  investigate  some  aspects  concerning  the 
illumination  of  an  artificial  earth  satellite  by  the  sun  in  relation  to 
the  track  of  the  sub- satellite  point  as  measured  in  a  geographic 
coordinate  system.     But  before  stating  the  complete  problem,   some 
preliminary  discussion  is  necessary. 

Consider  a  satellite  in  a  circular  orbit  around  the  earth.     Its 
position  may  be  completely  specified  at  any  time  (ignoring 
perturbations)  if  its  altitude,  orbital  inclination  to  the  equatorial 
plane,  and  the  longitude  of  its  ascending  node  at  a  particular  time  is 
known.     The  longitude  of  the  ascending  node  may  be  defined  as  the 
longitude  at  which  the  sub- satellite  point  crosses  the  equator  entering 
the  northern  hemisphere.     Under  these  circumstances  the  track  of 
the  sub- satellite  point  on  the  earth's  surface  may  be  computed  for  any 
set  of  sequential  orbits. 

Due  to  the  rotation  of  the  earth  about  its  polar  axis  ,  the  longitude 
of  the  ascending  node  is  displaced  westward  cars,  each  orbit  by  an 
amount  equal  to  the  rotational  velocity  of  the  earth  times  the  orbital 
period  of  the  satellite.     This  displacement  will  be  denoted  weT  ,  where 
we  is  the  angular  rotational  velocity  of  the  earth  about  its  polar  axis  , 
and  T  is  the  orbital  period  of  the  satellite. 

If  it  is  desired  to  determine  if  the  track  of  the  sub- satellite 
point  will  pass  within  a  given  distance  of  a  given  point  on  the  earth's 


surface  ,  then  under  the  stated  conditions  the  problem  is  completely- 
deterministic.     A  simple  but  somewhat  tedious  calculation  will 
ascertain  the  number  of  orbits  required,  the  exact  time,  and  the 
actual  distance  from  the  given  point  at  the  time  of  closest  approach. 

Now  suppose  that,  instead  of  being  completely  specified,  these 
parameters  of  altitude,  inclination,  and  longitude  of  ascending  node 
were  random  variables.     Conceivably  this  situation  could  arise  where 
a  satellite  has  been  "lost"   to  surface  detection  units  for  an  extended 
time  period.     Similarly,  this  proposition  might  be  required  in  a 
feasibility  study  of  a  ground  visual  detection  station.     In  any  event, 
under  these  circumstances  the  problem  of  ascertaining  if  the  satel- 
lite track  passes  within  a  given  distance  of  a  given  point  is  no  longer 
deterministic  but  stochastic.     The  solution  must  be  given  as  a 
probability  of  such  passage  after  a  given  number  of  orbits,  and  will 
be  a  function  of  these  three  random  variables. 

To  simplify  the  problem  considerably,  let  the  variables  of 
altitude  and  inclination  take  on  a  range  of  discrete  values.     Then  seek 
a  solution  for  the  random  variable  of  longitude  of  ascending  node  for 
each  pair  of  values  for  altitude  and  inclination.     The  longitudes  of 
the  second  and  all  subsequent  ascending  nodes  are  simple  functions 
of  the  first  nodal  longitude,   since  specifying  the  altitude  specifies 
orbital  period,  hence  specifies    the  longitude  displacement,  weT. 
Then  the  probability  of  passage  is  a  function  of  the  random  variable 


of  longitude  of  first  ascending  node,  and  of  the  orbital  period,   T. 

The  second  part  of  the  problem  is  to  determine  if  the  satellite 
is  also  illuminated  by  the  sun.     It  appears  obvious  that  the  satellite 
need  not  be  over  the  sunlit  side  of  the  earth  to  be  illuminated.     In 
fact,  the  only  requirement  for  illumination  is  that  the  satellite  be 
outside  the  shadow  cone  of  the  earth.     For  satellites  at  sufficient 
altitudes  this  may  not  occur  until  the  satellite  is  well  into  the  "back 
side"   of  its  orbit. 

With  these  considerations  in  mind,  the  main  problem  may  now 
be  stated: 

Given    a  satellite  at  an  altitude  h ,  with  an  orbital  inclination  i , 
and  two    geographical  points  on  the  earth's  surface  on  the  same 
meridian  of  longitude  W,  and  equidistant  from  the  equator  at  latitudes 
L     North  and  South: 

Find    the  joint  cumulative  probability  that  the  satellite  will  pass 
within  a  given  angular  distance  s    of  either  point  and  that  the  satellite 
will  simultaneously  be  illuminated  by  the  sun,  the  result  to  be  stated 
as  a  cumulative  probability  distribution  function  in  terms  of  successive 
orbital  periods  after  initial  reference  time  (arbitrarily  set  equal  to 
zero). 

Some  previous  studies  of  satellite  illumination  have  been  made, 
made,  but  no  studies  have  been  found  which  treat  the  problem  from  a 
stochastic  viewpoint. 


This  thesis  first  develops  the  necessary  mathematical  models 
for  the  orbital  mechanics  ,  then  sets  forth  the  precise  probability 
formulation  statements.     The  statistical  solutions  are  then  derived 
from  a  digital  computer  simulation  program  which  is  discussed  in 
the  appendices. 

While  the  design  and  writing  of  the  computer  simulation  program 
proved  of  great  interest  to  the  author,  and  provided  an  opportunity 
to  apply  recently  acquired  knowledge  of  Monte  Carlo  techniques,  it 
must  be  emphasized  that  the  computer  program  is  not  an  entity  in 
itself ,  but  merely  a  useful  tool  in  arriving  at  a  solution  to  the  stated 
problem. 


Leinf  elder  ,  E.J.  ,  and  Thrasher,  J.R.  ,  Orbit- Earth  Shadow 
Relations,  AST/EOR- 13277  ,  Chance-Vought  Corp.  ,  Dallas,  Texas, 
I960.       CONFIDENTIAL 

^Koelle  ,  H.H.  ,  Editor,  Handbook  of  Astronautical  Engineering, 
McGraw-Hill  Book  Co.  ,  New  York,    1961  . 


CHAPTER    I 
ORBITAL  MECHANICS 

Initial  Assumptions 

1.  The  earth's  solar  orbit  and  the  satellite  orbit  have  zero 
eccentricity. 

2.  The  earth  is  a  perfect  sphere  with  total  mass  concentrated 
at  the  center. 

3.  The  earth's  equatorial  plane  is  inclined  to  the  ecliptic 
plane  at  a  constant  angle  i'  (=  23°  26'    59") . 

4.  The  earth's  solar  umbra  is  a  circular  cylinder  (valid  for 
values  of  h  less  than  three  times  the  earth's  radius). 

5.  At  initial  reference  time  the  longitude  of  the  ascending  node 
of  the  satellite  may  be  considered  a  random  variable  ,  uniformly 
distributed  on  the  interval       0,2fT    .     At  this  time  the  satellite  will 

be  assumed  to  be  located  at  the  ascending  node. 

6.  All  other  perturbations  of  a  satellite  orbit  will  be  ignored. 
This  includes  atmospheric  drag  and  solar  radiation  pressure. 

Discussion  of  Assumptions 

The  effect  of  these  assumptions  is  to  confine  the  problem  to 
one  involving  only  central  forces.     This  also  implies  that  precession 
of  the  satellite  orbital  plane  will  not  be  included.     In  order  that  these 
omissions  do  not  radically  affect  the  solution,  the  probability  dis- 
tribution function  sought  will  only  be  examined  for  a  period  of  16  orbits 


This  will  be  on  the  order  of  24  hours  for  the  range  of  altitudes  con- 
sidered.    As  will  be  seen  subsequently,  two  more  random  variables 
must  be  introduced;  one  to  account  for  the  space  orientation  of  the 
satellite  orbital  plane,  and  another  to  account  for  seasonal  variations. 

Method  of  Attack 

A  mathematical  model  will  be  established  representing  the  motion 
of  the  satellite  and  the  sun  with  reference  to  a  geocentric  coordinate 
system,  followed  by  a  probability  formulation. 

Mathematical  Model 

Let  E    represent  the  vector  from  the  center  of  the  earth  to  the 
sun,  and  let  r    represent  the  vector  from  the  center  of  the  earth  to 
the  satellite.     The  initial  assumptions  will  ensure  that  these  vectors 
are  of  constant  magnitude  and  possess  constant  rotational  velocities 
about  the  center  of  the  earth. 

Using  Figure  I- 1 ,  find  the  direction  cosines  of  E     in  the  rectan- 
gular geocentric  coordinate  system  depicted,  where  the  X-axis  is 
along  the  vernal  equinox,  the  Z  axis  is  directed  through  the  north 
celestial  pole,  and  the  Y-axis  perpendicular  to  the  other  two  axes  in 
a  conventional  manner.     The  origin  is,  of  course,  the  center  of  the 
earth. 

During  a  yearly  cycle  the  angle  H  takes  on  all  values  in  the 
interval    10,  2ft     .     The  values  of  x,  y,  and  z  are  easily  determined  as: 


z 


Y 


Fig.  I-i:    Representation    of 
Earth-Sun    vector     E 


x  =  E  cos  H 

y  =  E  cos  i  '  sin  H 

z  =  E  sin  i  '  sin  H 
Hence  the  direction  cosines  are  :  (cos  H,  cos  i'  sin  H,   sin  i'  sinH). 

From  Figure  1-2  find  the  direction  cosines  of  r    in  the  same 
coordinate  system.     The  angle  D   may  be  thought  of  as  the  right 
ascension  of  the  ascending  node.     The  angle  G     represents  the 
angular  distance  of  r    from  the  ascending  node.     Obviously  the  angles 
D  and  G  take  on  all  value  in       0  ,  271    • 

From  Figure  1-2  the  angle  F   is  the  projection  of  G  onto  the 
equatorial  plane.     Let  r'  be  the  projection  of  r    onto  the  same  X-Y 
plane.     Then  r1  =  r  cos  (angle  between  r    and  r1). 

Using  Figure  1-3  the  right  spherical  triangle  with  sides  G,  F  and 
L   '  reveals  that  L  '  is  the  angle  between  r    and  r'.     Then 

sin  L'     =  sin  i  sin  G  (1-2) 

and 

cos  L'     «  (1  -  sin     L')  z 

=  (1  -  sin2  i     sin2    G)^  (1-3) 

Also  from  the  right  spherical  triangle  in  Figure  1-3: 

tan  F    =    cos  i  tan  G  (1-4) 

hence 

_  i 
F  =  tan       (cos  i  tan  G)  (1-5) 


z 


cx.sc ending   mode 


Fig.  1-2  :     Representation    of     E  ari 
Satellite     vector      F 


V 


F 

Fig.  I~3 :    Simplified     view    of 
Earth  -  Satellite     vector     r 
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Referring  again  to  Figure  1-2, 
x  =  r  '  cos  (D+F) 

=  r    cos  L  '  cos  (D+F) 

2  2         Vo 

=    r  (1  -  sin       i  sin     G)    z      cos 

y  =  r  "  sin  (D+F) 

/!     .  2  .   .  2  „.J4 

~  r  (1  -  sin       l  sin     G)  sin 

z  =  r  sin  L  ' 

=  r  sin  i  sin  G 


-1 

D  +  tan    (cos  i  tan  G) 

(1-6) 


D  +  tan"  (cos  i  tan  G) 
(1-7) 


(1-8) 


Using  Equations  (1-6),   (1-7)  and  (1-8),  the  direction  cosines  of  r    are: 


, ,     .  2  .  .  2  ^  x  v% 

(1  -  sm  i  sin     G)           cos 

/i         .    2  .  .    2  ^\^ 

(1  -  sm  i  sm     G)         sm 


D  +  tan"    (cos  i  tan  G) 
-1 


D  +  tan       (cos  i  tan  G) 

(1-9) 

Referring  to  Figure  1-4  it  is  noted  that  the  angle  G  describes 
the  angular  position  of  the  satellite  (or  the  sub- satellite  point)  in  the 
orbital  plane  from  the  ascending  node  at  time  t.     Hence 


G  -  w_t 

s 


1-10) 


where 


w        =  angular  velocity  of  satellite  in  radians/unit  time 
t  -  instantaneous  time  after  reference  time. 

As  the  satellite  proceeds  in  its  orbit  to  time  t,  the  earth  also  rotates 
about  its  axis  with  an  angular  velocity  we  .  Then  dl  is  the  apparent 
change  in  the  longitude  of  the  sub-satellite  point  in  time  t.     Clearly, 
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Fig.  I-^t:    Location    of    sub~  satellite 
point     at    time     t 
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F    =  wet  +  dl         .  (1-11) 

Then  the  longitude  of  the  sub- satellite  point  at  time  t  is: 

l't     =  l'o    +dl 

=     l'o    +  F  -  wet 

=     1'0    +  tan-1  (cos  i  tan  G)     —  wgt  (1-12) 

It  is  more  convenient  at  this  point  to  develop  an  analytic  model 
which  will  indicate  whether  or  not  the  satellite  is  illuminated  by  the 
sun.     Figure  1-5  depicts  the  plane  containing  the  E    and  r    vectors. 
For  simplicity,  let  both  E    and  r    be  unit  vectors. 

As  the  satellite  orbits  about  the  earth,  the  angle  A     (between 
the  E    and  r    vectors)  can  take  on  all  values  in  the  interval       0  ,Zfx\  . 
However,   since  the  representation  in  Figure  1-5  is  symmetric  about 
the  E    vector,  the  problem  requires  an  interest  in  the  magnitude  of 
A  only  in  the  interval       0  ,  7T  J    . 
Clearly , 

E  •     7      =    cos  A  (1-13) 

If  cos  A  is  non-negative,  the  satellite  will  be  illuminated  by  the  sun. 
If  cos  A  is  negative,  the  satellite  may  or  may  not  be  illuminated. 
Let  B      =    A  -    Tt  /2  as  shown  in  Figure  1-5.     Let  C    be  the  angle  at 
the  satellite  subtended  by  the  center  of  the  earth  and  the  surface  of  the 
earth.     Then 

C  -  sin"1    /  Re 

Re    +h 
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satellite 


Sun 


Fig. J? 5:    Plane    containing   E  and  T  vectors 


14 


a  constant  for  circular  orbits,   where  Re    is  the  radius  of  the  earth. 

Then  if  B  +  C  is  greater  than      7T/2,  the  satellite  will  be  in  the 

shadow  cylinder  of  the  earth.     Hence  if  the  satellite  is  illuminated, 

then 

cos        (E  •  7)  7t       +  sin"       /  Re      )  (1-15) 

Re      +h 


or 


i     __  1      /  R  \      / 

cos      (E     •     r)    +  sin"        e        1    =    TT      .  (1-16) 


Re    +h 


If  it  is  not  illuminated,  then 


cos_1(E  •     r)    +  sin"1    [  Re  [     >    j(  (1-17) 

\Re      +h       / 

Using  the  results  of  Equations  (1-1)  and  1-9) ,  then: 

,      2  7  '/ 

cos  A  -  E  •   r  =  (1  -  sin      i  sin      G)  x 

<  cos  H  cos        D  +  tan"    (cos  i  tan  G)         + 

cos  i  '     sin  H  sin      D  +  tan"    (cos  i  tan  G)      > 

+  sin  i  '     sin  i  sin  H  sin  G  (1-18) 

The  next  stage  of  the  problem  is  to  determine  if  the  satellite 
track  passes  within  a  given  angular  distance  s  of  either  of  two 
geographical  points  with  coordinates  (L,   W)  and  (-L,,  W).     With  this 
concept  the  satellite  track  may  be  considered  to  sweep  out  a  "swath" 
over  the  earth's  surface,  of  width  2s.     Denote   1 '0    as  the  longitude 
of  the  ascending  node. 
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For  illustrative  purposes  it  is  advantageous  to  show  a  schematic 
plot  of  the  satellite  track  across  the  surface  of  the  earth,  with 
special  emphasis  upon  its  relationship  to  the  designated  geographical 
points.     For  clarity,  a  Mercator  projection  will  be  used. 


For  cases  where 


-  s     > 


then  the  situation  in 


Figure  1-6  will  exist.     Clearly,  if  the  value  of  W  lies  in  one  of  the 
intervals  lj  ,     h  +  i\      for  j  =   1,   3,   5,   7,  then  the  desired  passage 


will  occur.     Indeed  this  is  a  necessary  condition. 


For  the  case  where 


-   s 


+  s  ,  then  Figure 


1-7  will  apply.     Here,  the  desired  passage  will  occur  if  W  lies  in 
one  of  the  intervals       1;  ,     li-f  i         >     for  j     =     1,3.     This  too  is  a 
necessary  condition. 

Using  these  conditions  to  determine  if  passage  occurs,  the  time 
t  of  closest  approach  to  the  geographical  point  may  be  evaluated. 
Using  this  value  of  t  to  compute  cos  A  in  Equation  (1-18),  then  the 
conditions  contained  in  Equation  (1-16)  are  applied  to  determine  if 
the  desired  illumination  occurs. 

It  is  now  necessary  to  develop  expressions  for  the  orbital  period 

and  angular  velocities  for  the  satellite  and  the  earth.     It  can  be  shown 

that  the  period  of  a  satellite  is: 

3/2 


2  7t    (a); 

(g'm')1/2 


.1-19) 
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where 


a        =  semi-major  axis  of  the  orbital  ellipse 

=  Re  +  h,  for  circular  orbits 

/ 
G        =  universal  gravitational  constant 

M       =  mass  of  the  earth. 


In  mks  units  , 


Re     =  6.371     X    10°    meters 

G        =  6.668     X  10"  newton-meter    /kilogram 


Then 


'  24 

M       =5.98    X    10        kilograms 


2    71   (6.371     X     10b  meters         +  h)~" 

6.668  X  10'11     nt.  -mZ    X5.98X1024kg 
kg2 

=     2  71(6.371     X     106  m.     +  h) 


1/2 


(3.99    X     10 M    m3  /sec2)  1/2 
=  3.  148  X  10~7  (6.371  X  106    +  h)  3/2    sec 
-5.246    X  (6.371     +  h)3/2    min.  (1-20) 


where  h  is  given  in  units  of  10       meters. 
Since  w  aT     =  2  Tt       , 


and 


wg    =     2  tt    /T  (1-21) 


we    =  2  n   /Te 


where  T       is  the  rotational  period  of  the  earth  and 
e  r 
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we     =    (2  n  /1440)  min. 

=    4.363    X    10~3  min.  (1-22) 
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CHAPTER  II 
PROBABILITY    FORMULATION 

The  next  stage  leading  to  a  solution  of  this  problem  is  the 
formulation  of  a  probability  statement.     The  probability  sought  is 
the  probability  that  illumination  occurs  simultaneously  with  the 
desired  passage  by  the  n        orbit,  for  n  =  1,   2  ,    .  .  .  ,    16. 

Looking  at  the  first  orbit,  the  longitude  system  may,  without  any 
loss  of  generality,  be  arbitrarily  defined  such  that  the  longitude  of  the 
ascending  node,    1'0   ,  is  equal  to  zero,  and  all  other  longitudes  are 
measured  eastward  from  zero  through  a  total  angle  of  2  7T   . 

Let  M  be  the  event  that  the  desired  passage  occurs  sometime 
during  the  first  orbit,  and  let  N  be  the  event  that  the  desired 
illumination  occurs  at  the  particular  time  t  when  the  satellite  makes 
its  closest  approach  to  the  geographical  point.     Note  that  M  is  a 
time-independent  event,  while  N  is  strictly  dependent  upon  the  value 
of  t.     The  probability  sought  is  P   I  MN 

Studying  each  event  separately,  it  is  clear  that 

P    [  m]      =  P  [  lj      =  W      =  L  +  1]  (2-1) 

where  j  =  1  and  3,  or  else  j  =  1  ,   3,   5,7,  depending  upon  the  values 
of  i  ,   s  and  L. 

This  equation  can  be  written  in  more  meaningful  language  if 
P       Mi      is  defined  as  the  probability  that  W  lies  in        L,     lj  +  lj 
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Clearly,  the  events  Mj,  M3  ,    .  .  .  ,  are  mutually  exclusive, 
hence 

P 


M  -  Z  p  tMj] 


(2-2) 


odd  j 


Also,  from  Equation  (1-16), 

p[n]-p[a+C    =ftl     =  P    [  A    =  It   -     C]  (2-3) 

Hence  the  random  variable  associated  with  event  M  is  W,  and 

it  appears  intuitively  correct  that  A  is  the  random  variable  associated 

with  event  N.     Then  the  joint  probability  sought  is 

P    f  MN]      =         y  P  MN  (2-4) 

odd  j 

To  verify  that  A  is  a  random  variable,  note  that  Equation  (1-18) 
shows  that  A  is  a  function  of  i  ,  G,  D  and  H.     But  for  a  particular 
time  t,  the  values  of  i  and  G  are  constant.     Then  PA      =  71  -  C 
is  simply  the  probability  that  D  and  H  take  on  a  range  of  values  such 
that  A     =  7T  —    C  ,  which  implies  that  D  and  H  are  random  variables 
with  associated  probability  density  functions,  fj^(d)  and  fjj(h) ,  both 


defined  on  the  interval        0  ,   2  7T 


[0, 


Care  should  be  taken  to  avoid 
confusion  with  h,  the  satellite  altitude;  however,  in  most  situations 
the  meaning  will  be  clear  from  the  context. 

Using  the  basic  assumptions  and  the  physical  interpretations  of 
D  and  H,  it  appears  intuitively  correct  to  state  that  they  are  inde- 
pendent, uniformly  distributed  random  variables  such  that 
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fD(d)  1/271 


fH(h) 


0 


<         < 

for  0    =  d    -  2tz 


elsewhere 


for  0    =  h    *    2n\ 


elsewhere 


(2-5) 


(2-6) 


It  quickly  follows  that  A  =  A(G,H)  is  a  randorn  variable  with 
an  associated  probability  density  function,  f«(a),  defined  over  the 
interval       0  ,  7T     . 

While  it  appears  unlikely  that  the  mathematical  expression  for 
f  a  (a)  could  be  derived  quickly  by  analytic  methods  ,  the  formulation 
of  the  problem  will  continue  as  if  f »  (a)  were  completely  described. 

Before  proceeding  further  with  Equation  (2-4)  it  is  necessary  to 

develop  the  form  of  the  probability  density  function  for  W.     Since  the 

only  quantity  of  interest  is  W  -   1'q    ,  and  since  one  of  the  basic  initial 

assumptions  was  that  1'0    was  uniformly  distributed,  it  follows  that 

W  is  also  uniformly  distributed.     Arbitrarily  setting  l'G    equal  to 

zero,  as  before,  merely  ensures  that  W  is  defined  over  the  interval 

[   0,2  71  ]. 

Then 

/ 

fw(w)    =  1/2 -ft  for  0     =  w  =  2t 

0  elsewhere 


(2-7) 


Now  by  the  rules  for  conditional  probabilities, 


M.N 
J 


=  P 


N 


M- 


M; 


(2-8) 
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and 


MN 


N 


M 


oddj 


J 


M 


J 


(2-9) 


Assume  now  that ,  associated  with  P    N  M-       ,  there  exists  a 
conditional  probability  density  function  f  *  |w(a    w) »  defined  over  the 
,  K      .       Then  Equation  (2-9)  may  be  written  as 


interval 


0 


MN      = 


odd  i       1-  0 

J        J 


"AW 


[aw)  f™(w)  da    dw 


1 
2  7C 


^j  +  1        '^ 

I  /   / 


lj+1      *c 

f^|w(a|w)    da  dw        (2-10) 
oddj        L  0 

using  the  results  of  Equation  (2-7). 

The  final  probability  to  be  calculated  is  the  joint  cumulative 
probability  that  the  desired  events  M  and  N  will  occur  simultaneously 

at  least  once  by  the  n        orbit,  for  n  =  1,  2 16.     To  accomplish 

this  requires  some  new  notation. 

Let  P  represent  the  desired  probability  of  at  least  one  such 
joint  occurrence  by  the  n  orbit.  Let  P(l)  be  the  probability  that 
the  desired  event  occurs  on  the  first  orbit;  clearly  P(i)  is  the  solution 

to  Equation  (2-1.0).     Let  P(n' ,n+l)  be  the  joint  probability  of  no  such 

st 

occurrence  on  the  first  n  orbits  and  one  occurrence  on  the  (nil) 

orbit. 
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Then 


Px      =  P(l) 


P2     -Pj  +  Pfl'^) 
P3     -P2  +  P(2,,3) 
and,  in  general, 

Pn+1    =  Pn    +    F(n',n+1)       •  (2-H) 

This  equation  describes  symbolically  the  procedure  for 
calculating  Pn    for  n  =  1,  2,    .  .  .  ,     16.     The  problem  of  calculating 
the  probabilities  for  orbits  subsequent  to  the  first  can  be  considerably 
alleviated  by  utilizing  a  characteristic  of  the  satellite  tracks. 

From  Equation  (1-12)  it  is  easily  shown  that  the  longitude  of  the 
ascending  node  for  the  second  orbit  is 

rT     =     Vo        "WeT       '  <2-12> 

i.e.  ,    1'iji    has  been  displaced  westward  of  1'       by  the  amount  w    T. 
But  this  is  equivalent  to  retaining   1'0    as  the  longitude  of  the  ascending 
node  for  subsequent  orbits  and  displacing  W  eastward  a  distance  w    T. 
In  other  words,   Equation  (2-10)  may  be  used  for  the  n^n    orbit,  pro- 
viding the  longitude  of  the  geographical  points  is  set  equal  to 

(W)n    =■:  W  +  (n  -  1)  wgT      for  n  =   1,2,    .  .  .  ,    16  (2-13) 

Nonetheless,  there  remains  unsolved  the  problems  of  determining 
the  limits  of  integration  (with  respect  to  w)  and  discovering  the 
analytic  expression  for  f .  |w(a|w)  in  Equation  (2-10).     These  problems 
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will  be  avoided  by  abandoning  the  analytic  approach  in  ia-,-c:t 
of  a  computer  simulation  program  using  Monte  Carlo  methods 


26 


CHAPTER  III 
SIMULATION  PROBLEM 

It  is  expedient  first  to  examine  the  form  of  the  desired  solution 
in  order  to  assess  the  proper  method  of  arriving  at  it. 

The  cumulative  probability  distribution  function  sought  is  not  a 
continuous  function,  but  a  monotonic  step  functions  with  certain 
probabilities  associated  with  each  point  corresponding  to  the  end  of 
an  orbital  period.     In  fact,  the  desired  result  will  be  a  distribution 
function  truncated  at  16  orbital  periods. 

Now  assume  for  the  moment  that  the  problem  can  be  set  up  for 
simulation,  where  a  suitable  number,  k,  of  samples  can  be  obtained 
from  a  pseudo-random  number  generator.     Let 

Xi    =  number  of  times  that  passage  occurs  on  the  1st  orbit 

Xj     =  number  of  times  that  passage  occurs  on  the  itn    orbit 

and  not  on  the  i— 1  previous  orbits  ,  for  i  =  2  ,3  , ...  16. 

st 
y-i    =  number  of  times  that  illumination  occurs  on  the   1        orbit, 

given  that  simultaneous  passage  occurs  on  the  first  orbit. 

y^    =   number  of  times  that  illumination  occurs  on  the  i        orbit, 

given  simultaneous  passage  on  the  same  orbit,  and  that 

such  joint  event  did  not  occur  on  the  i— 1  previous  orbits  , 

f  o  r  i  =  2 , 3 , . . .  16 . 

Here  ,   "passage"   implies  that  the  satellite  track  passes  within  the 
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specified  angular  distance  s  of  either  of  the  two  geographical  points 
With  these  precise  definitions  it  is  possible  to  assert;  that: 


1 


E 


['*] 


=      E 


y- 


p(i) 


yi 


+       E 


[P(l',2)] 


yn    +  y« 


and,  in  general , 


1 


Zn 
i-1 


yi 


(3-1) 


If  P   (p)  is  defined  to  be  the  cumulative  probability  of  passage 


(independent  of  illumination),  then,  in  similar  manner, 
E     [Pn(p)]      ~     iV  x. 

kZ_,i=1 


Now  define  u-;    and  v--     such  that 


(3-2) 


u 


1J 


V 


1J 


'  1        if  the  desired  passage  occurs  for      \ 


the  first  time  on  the  i  x     orbit, 

tVi 
for  the  j        sample  random  number 

value 

^  0        Otherwise 

(    1        If  the  desired  illumination  and 

passage  occurs  for  the  first  time 
on  the  i       orbit ,  for  the  -y—  sam- 
ple random  number  value. 
3        Otherwise 


M3-3) 


>  (3-4) 
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where  j  =  1,2,.  .  .  ,k.     Note  that  v--     =  1  does  not  require  that  u^;     =  1. 


Then  k 


k 


(3-5) 


)' 

)■ 

V.  . 

1J 

(3-7) 

/ 

1=1 

L i 

j  =  l 

n 
v ' 

k 

1 

k 

7, 

y 

uij 

(3-8) 

and  y.  =        £_,  v^  (3-6) 


Hence  n 

and 

E     [Pn(p>]       = 

K. 

i=l  j=l 

Equations  (3-7)  and  (3-8)  can  be  considered  as  the  basis  for  the 
simulation  model. 

The  design  of  the  computer  simulation  program  is  described  in 
some  detail  in  Appendix  A  ,  with  a  description  of  some  of  the  pro- 
cedures used  to  reduce  the  required  computation  time.     The  language 
used  in  writing  this  program  was  the  CDC-1604  version  of  NELIAC  , 
which  is  peculiarly  well-adapted  to  this  type  of  program. 

Appendix  B  contains  a  summary  of  the  output  of  this  program, 

along  with  a  sensitivity  analysis  for  the  effect  of  single  parametric 

variations  of  satellite  altitude,  orbital  inclination,  latitude  of  geo- 
graphical points,  and  angular  distance  s. 

29 


BIBLIOGRAPHY 


1.  Space  Navigation  Handbook,  NAVPERS  92988. 

2.  Merrill,  G.  ,   Editor,     Principles  of  Guided  Missile  Design, 
D.   Van  Nostrand  Co.  ,  Inc.  ,  Princeton,  New  Jersey,    I960. 

3.  Van  Allen,  J. A.  ,  Scientific  Uses  of  Earth  Satellites, 
University  of  Michigan  Press,  Ann  Arbor,  Mich.  ,  1956. 

4.  The  NELIAC  Compiler  Language,  Departmep.t  of  Operations 
Research,   U.S.   Naval  Postgraduate  School,  Monterey,  Calif.  , 
1963. 


30 


APPENDIX    A 
DESIGN  OF  A  COMPUTER  SIMULATION  PROGRAM 

The  primary  concept  of  the  program  is  an  artificial-time - 
controlled  sequence  of  computations,  based  upon  Equation  (1-16)  and 
a  modification  of  Figures  1-6  and  1-7.     The  satellite  is  figuratively 
put  in  orbit  at  its  ascending  node,   corresponding  to  time  t  =  0,  and 
where  the  longitude   1'       is  also  identically  zero.     A  suitable  time 
increment  is  used  to  activate  an  artificial  time  "clock"    which 
determines  the  instantaneous  latitude  and  longitude  of  the  sub- 
satellite  point. 

Certain  short-cuts  are  immediately  available  to  reduce  the 
otherwise  inordinate  number  of  calculations  required. 

Figure  A- lis  similar  to  Figure  1-6,   except  that  the  end  points  of 
the  longitude  intervals  correspond  to  the  instantaneous  longitudes  of 
the  satellite  track  when  the  instantaneous  latitudes  equal  either  L  -  s 
or  L  +  s .     Clearly  the  desired  passage  cannot  occur  unless  the 
instantaneous  track  latitude  L't    lies  in  the  interval    [L-s,   L+sJ        , 
which  implies  that  if  the  instantaneous  longitude   l't    lies  in  the  interval 

L,   liAj     for  j  =  1,  3  ,  5 ,  7 ,  then  passage  will  occur  if  W  also  lies 
within  this  same  interval. 

In  similar  manner,  the  interval     I  L  ,    L+3J    for  j    -  1,   5,  is  defined 
for  situations  where  L  -  s  ~  i    =    L+s.     See  Figure  A -2.     Again, 
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clearly,  passage  will  occur  if  l't    and  W  lie  within  this  interval.     The 

numbering  of  the  subscript  j  is  designed  to  facilitate  the  computer 

programming . 

From  Equations  (1-2)  and  (1-10): 

sin  L'     =     sin  i  sin  w0t 

hence  for  L1    =    L-s, 

t    =     1         sin"1    |      sin  (L   — s) 
w  I  sin  i 

and  for  L't    =  L+s  , 

t    =   1  sin"1    J   sin  (L+s)  ]  -t2  (A-2) 

ws  ^         sin  i 

where  t:     corresponds  to  the  time  at  which  the  satellite  track  attains 
a  value  of  longitude  equal  to  h.     The  remaining  values  of  t-  for  j 
2  can  be  determined  quickly  by 

t3      =        T       -        t2 

2 

*4   =    i   -    h 

2 
t5  T       +       tj 

2 

*6      =       1      +       h 

2 

t?       =        T        +        t3 

2 
tg       =        T        +        t4 

2 

which  can  be  easily  verified  by  the  symmetry  of  the  track.     T  is  of 
course  the  satellite  period. 

Using  Equation  (1-12)  it  follows  that 

A-4 


1.     =  tan     (cos  i  tan  w„ti)     -w„t. 
J  s  j'  e  j 


and  accordingly 

lj    =  tan_i(cos  i  tan  wgt])     -w  tj 
12    =  tan_1(cos  i  tan  wst2j     .  w  t2 


(A -3) 

(A  -4) 


where   lj  and  L,    are  defined  as  in  Figures  A-land  A-2.     In  these 
figures  if  W  lies  within  any  of  the  intervals  depicted  then  passage 
will  occur.     The  remaining  values  of  1-     can  be  computed  in  exactly 
the  same  manner  as  Equations  (A -3)  and  (A-l). 

Then  the  first  test  for  the  desired  passage  is  determining  if  the 
sample  value  for  the  random  variable  W  lies  within  one  of  the  afore- 
mentioned intervals.     If  W  satisfies  this  test,  then  the  test  for  illumi- 
nation may  be  made.     However,  it  is  possible  for  W  to  be  outside 
such  intervals  and  still  provide  for  the  desired  passage.     If  W  lies 
outside  of,  but  sufficiently  close  to,   either  end  point  of  the  depicted 
intervals,  then  passage  may  occur  as  in  Figure  A-3. 

At  this  point  it  is  appropriate  to  introduce  the  approximation 
s  '  of  the  true  instantaneous  angular  distance  between  the  sub- satellite 
point  and  the  designated  geographical  point.     For  values  of  s1 
sufficiently  small,  the  curved  surface  of  the  earth  can  be  approximated 
by  a  plane.     Then,  from  Figure  A-3, 


s'  = 


IA 


W    -  1' 


2  2t  . 

cos    L\ 


1/2 


(A -5) 


A-5 


W       1;        i\ 


"Satellite 
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Fi.g.A'3: 


ercator    projection    of    trac«; 
special      case 
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where  the  cosine  term  is  used  to  compensate  for  the  "shrinkage" 
in  distance  between  adjacent  meridians  of  longitude  with  increasing 
latitude . 

(w  -rt)  x 


=  s  ,  hence  this  quantity  may  be  used  as  an  upper  bound 


It  should  be  obvious  that  passage  cannot  occur  if 
cos  L't 

in  an  iterative  procedure  for  computing  successive  values  of  s'  for 
successive  values  of  time.     Due  to  the  monotonic  nature  of  the 
satellite  track,  a  minimum  value  of  s'  will  exist  for  some  value  of  t, 
say  t*.     This  value  of  s"  is  compared  to  s;  if  s1     =  s,  then  the  desired 
passage  will  occur.     Otherwise,  passage  cannot  occur. 

If  these  tests  are  applied  to  the  i        orbit  and  j       sample  value  of 
the  random  variable  W,  the  value  of  ujj    in  Equation  (3-3)  is  determined 
If  uj:     =  1,  then  the  illumination  condition  may  be  determined  by 

J 

application  of  Equation  (1-16).     Applying  this  test  to  the  i       orbit  and 
j         sample  values  of  the  random  variables  W,  D  and  H  determines 
the  value  of  vjj    in  Equation  (3-4). 

However,  it  may  happen  that  the  desired  passage  will  occur  with 
no  illumination  on  an  earlier  orbit,  but  with  the  desired  passage  and 
illumination  on  a  later  orbit.     In  order  to  compute  Pn(p) ,  it  is 
necessary  to  record  the  earlier  event  where  U£j     =     1,  but  continuing 
the  process  in  order  to  compute  Pn. 

Since  this  simulation  program  was  only  incidental  to  the  solution 
of  the  main  problem,  the  computer  flow  diagrams  and  NELIAC 
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language  flow  charts  have  not  been  included.     However,  the  IBM 
binary  card  decks  and  program  print-out  are  available  from 
Professor  Richard  M.   Thatcher,  of  the  Operations  Analysis 
Department. 
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APPENDIX  B 
SUMMARY  OF  PROGRAM  OUTPUT 

The  parametric  inputs  to  the  simulation  program  were  satellite 
altitude  h,  orbital  inclination  i,  latitude  of  designated  geographical 
points  L,  and  angular  distance  from  these  points  s.     The  values 
of  h  and  s  are  given  in  units  of  nautical  miles  ,  but  were  converted 
internally  to  units  of  10       meters  and  radians  respectively;  the 
values  of  i  and  L  are  given  in  degrees  ,  but  were  likewise  converted 
internally  to  radians.     The  values  of  these  parameters  and  the  in- 
cremental values  used  in  the  program  are  given  below: 

h:  100(300)1900  nautical  miles 

i:  20°(8°)68o 

L:  25°(10°)65° 

s:  100(100)500  nautical  miles 

The  value  of  k  (sample  size)  was  hopefully  selected  to  be  200  , 
in  order  that  the  approximations  contained  in  Chapter  III  might  be 
statistically  valid,  and  yet  require  a  reasonably  short  amount  of 
computer  time.     Even  so,  about  three  and  one-half  hours  actual 
operating  time  was  consumed  on  the  computer  for  a  single  run 
through  all  values  of  the  above  parameters. 

One  immediate  result  from  the  output  of  this  program  was  the 
discovery  of  a  "semi-optimal"  distance  s  =  300  nautical  miles. 
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In  many  instances,  for  various  values  of  h,  i  and  L,  the  cumulative 
probabilities  revealed  little  or  no  increase  with  values  of  s  greater 
than  300.     In  a  significant  number  of  cases  the  rate  of  increase  in 
probability  declined  sharply  after  s  =  300;    it  is  to  be  noted,  however, 
that  this  was  not  always  the  case.     Frequently  the  rate  of  increase 
in  probability  with  increasing  s  was  approximately  linear. 
Nevertheless  ,  it  was  decided  to  limit  the  tabulated  values  of  the 
probabilities  to  those  computed  for  s  =  300.     To  have  done  otherwise 
would  have  increased  the  size  of  this  appendix  by  a  factor  of  five. 

Another  unexpected  discovery  concerned  certain  substantial 
decreases  in  cumulative  probabilities  for  h  =  1900  nautical  miles. 
The  reduction  did  not  occur  consistently  for  all  values  of  the  other 
parameters,  however.     It  is  suspected  that  this  reduction  is  most 
likely  a  statistical  anomaly;  unfortunately  insufficient  time  was 
available  for  further  repeated  runs  using  this  value  of  h.     This  so- 
called  anomaly  is  suggested  as  a  topic  for  further  study,  as  is  an 
analysis  of  variance  for  repeated  runs.     The  computer  program  is 
easily  adapted  to  the  study  of  the  effects  of  variation  of  any  one 
parameter. 

Figures  B-l  through  B-10  give  the  joint  probability  of  passage 
and  illumination  vs.  number  of  orbital  periods,  while  Figures 
B- 11  through  B-15  give  the  independent  cumulative  probability  of 
passage  vs.   number  of  orbital  periods. 
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Figures  B-l  through  B-4  and  B-ll  through  B-15  represent 
probabilities  averaged  for  all  the  altitudes  examined  and,  as  such, 
show  the  sensitivity  due  to  variations  in  inclination.     The  rationale 
for  this  averaging  process  lies  in  the  fact  that  the  variations  in 
probabilities  due  to  increasing  altitude  were  sufficiently  small  to 
warrant  doing  so. 

Figures  B-5  through  B-10  ,  however,  are  given  for  specific 
latitudes  and  inclinations,  and  show  the  sensitivities  due  to 
variations  in  altitude.     In  these  cases  it  will  be  noted  that  the  values 
of  L  and  i  are  approximately  the  same.     In  these  cases  the  probabi- 
lities are  consistently  higher  than  in   cases  where  i  is  considerably 
greater  than  L.     Of  course,  when  i  is  less  than  L  by  3       or  more, 
no  passage  occurs. 

These  graphs  have  been  smoothed  to  some  extent  over  the  raw 
statistical  output,  but  the  general  trend  has  been  retained  in  all 
cases.     No  attempt  should  be  made  to  interpolate  between  discrete 
orbital  periods. 
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